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Abstract. The goal of this paper is to extend independent subspace 
analysis (ISA) to the case of (i) nonparametric, not strictly station- 
ary source dynamics and (ii) unknown source component dimensions. 
We make use of functional autoregressive (fAR) processes to model the 
temporal evolution of the hidden sources. An extension of the ISA sep- 
aration principle-which states that the ISA problem can be solved by 
traditional independent component analysis (ICA) and clustering of the 
ICA elements-is derived for the solution of the defined fAR independent 
process analysis task (fAR-IPA): applying fAR identification we reduce 
the problem to ISA. A local averaging approach, the Nadaraya- Watson 
kernel regression technique is adapted to obtain strongly consistent fAR 
estimation. We extend the Amari-index to different dimensional compo- 
nents and illustrate the efficiency of the fAR-IPA approach by numerical 
examples. 

Key words: nonparametric source dynamics, separation principle, ker- 
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1 Introduction 

Independent Component Analysis (ICA) [1,2] has received considerable attention 
in signal processing. One may consider ICA as a cocktail party problem: we 
have D speakers (sources) and D microphones (sensors), which measure the 
mixed signals emitted by the sources. The task is to recover the original sources 
from the mixed observations only. For a recent review about ICA, see [3]. In 
ICA the hidden independent sources are one-dimensional. The model is more 
realistic if we assume that not all, but only some groups of the hidden sources are 
independent ('speakers are talking in groups'). This is the Independent Subspace 
Analysis (ISA) generalization of the ICA problem [4] . The ISA model has already 
had some exciting applications including (i) the analysis of EEC, fMRI, ECG 
signals and gene data, (ii) pattern and face direction recognition. For a recent 
review of ISA techniques, see [5]. 
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One can relax the traditional independent identically distributed (i.i.d.) as- 
sumption of ISA and model the temporal evolution of the sources, for example, 
by autoregressive (AR) processes [6], however to the best of our knowledge, the 
general case of sources with unknown, nonparametric dynamics has been hardly 
touched in the literature [7,8]. [8] focused on the separation of stationary and er- 
godic source components of known and equal dimensions, in case of constrained 
mixing matrices. [7] was dealing with wide sense stationary sources that (i) are 
supposed to be block-decorrelated for all time-shifts, and (ii) have equal and 
known dimensional source components. 

One of the most exciting and fundamental hypotheses of the ICA research 
is due to Jean-Francois Cardoso, who conjectured that the solution of the ISA 
problem can be separated [4] to (i) applying traditional ICA and then (ii) clus- 
tering of the ICA elements into statistically dependent groups. While the extent 
of this conjecture, the ISA separation principle is still on open issue, it has been 
rigorously proven for some distribution types [9]. The goal of the present pa- 
per is to address the problem of ISA with nonparametric dynamics. Beyond the 
extension to not necessarily stationary dynamics, the temporal evolution of the 
sources can be coupled (it is sufficient that their driving noises are independent, 
no block-decorrelatedness for the sources are required) and we treat the case of 
unknown source component dimensions. We model the dynamics of the sources 
by functional AR (fAR) processes and derive separation principle based solu- 
tion for the resulting problem: the task is transformed to fAR estimation and 
ISA. To obtain strongly consistent fAR estimation the Nadaraya- Watson kernel 
regression technique is invoked. 

The paper is structured as follows: Section 2 formulates the problem domain. 
Section 3 shows how to transform the problem to functional AR estimation task 
and ISA, and presents the kernel regression based approach. Section 4 contains 
the numerical illustrations. Conclusions are drawn in Section 5. 



2 The Functional Autoregressive Independent Process 
Analysis Model 

We define the functional autoregressive independent process analysis (fAR-IPA) 
model. Let us assume that the observation (x) is linear mixture (A) of the 
hidden source (s), which evolves according to an unknown fAR dynamics (f) 
with independent driving noises (e). Formally, 

s t = f(s t _i, . . . ,s t - p ) +e t , (1) 
x t =As t , (2) 

where the unknown mixing matrix A 6 R DxD is invertible, p is the order 
of the process and the e m £ M. dm components of e = [e 1 ; . . . ; e M ] € M. D 

(D = ^2 m= idm) are (i) non-Gaussian, (ii) i.i.d. in time t and (iii) indepen- 
dent, I(e 1 , . . . ,e M ) = 0, where I denotes mutual information. The goal of the 
fAR-IPA problem is to estimate (i) the inverse of mixing matrix A, W = A -1 
and (ii) the original source s t by using observations x t only. 
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3 Method 

The estimation of the fAR-IPA problem (l)-(2) can be accomplished as follows. 
The observation process x is invertiblc linear transformation of the hidden fAR 
source process s and thus it is also fAR process with innovation Ae t 

x t = As t = Af (s f _i, . . . , s t _p) + Ae t = (3) 
= Af {A-^t-i, • ■ ■ , A _1 x t _ p ) + Ae t = (4) 
= g(x t _i,...,Xi_p) +n t , (5) 

where function g(ui, . . . , u p ) = Af (A _1 ui, . . . , A _1 u p ) describes the temporal 
evolution of x and n t = Ae t stands for the driving noise of the observation. 
Making use of this form, the fAR-IPA estimation can be carried out by fAR fit 
to observation x followed by ISA on n t , the estimated innovation of x. 

Let us notice that Eq. (5) can be considered as a nonparametric regression 
problem, we have u t = [x t _i, . . . , x t _ p ], v t = x f (t = 1, . . . ,T) samples from the 
unknown relation 

v t = g(u t ) + n 4 , (6) 

where u, v and n is the explanatory-, response variable and noise, respec- 
tively, and g is the unknown conditional mean or regression function. Non- 
parametric techniques can be applied to estimate the unknown mean function 
g(U) = i^(V|U), e.g., by carrying out kernel density estimation for random 
variables (u,v) and u, where E stands for expectation. The resulting Nadaraya- 
Watson estimator (i) takes the simple form 

where K and h > denotes the applied kernel (a non-negative real- valued func- 
tion that integrates to one) and bandwith, respectively, and (ii) can be used to 
provide a strongly consistent estimation of the regression function g for station- 
ary x processes [10]. It has been shown recently [11], that for first order and 
only asymptotically stationary fAR processes, under mild regularity conditions, 
one can get strongly constistent estimation for innovation n by applying the 
recursive version of the Nadaraya- Watson estimator 

«„) . EL"W(u-u,)) 

Er.,f D K(i"(u-u,)) 

where the bandwith is parameterized by (3 € (0, 1/-D). 



4 Illustrations 

Now we illustrate the efficiency of the algorithm presented in Section 3. Test 
databases are described in Section 4.1. To evaluate the solutions we use the 
performance measure given in Section 4.2. The numerical results are summarized 
in Section 4.3. 
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4.1 Databases 

We define three databases to study our identification algorithm. The smiley test 
has 2-dimensional source components [d m = 2) generated from images of the 6 
basic facial expressions 1 , see Fig. 1(a). Sources e m were generated by sampling 
2-dimensional coordinates proportional to the corresponding pixel intensities. 
In other words, the 2-dimensional images were considered as density functions. 
M < 6 was chosen. In the d-geom dataset e m s were random variables uniformly 
distributed on <i m -dimensioiial geometric forms. Geometrical forms were chosen 
as follows. We used: (i) the surface of the unit ball, (ii) the straight lines that 
connect the opposing corners of the unit cube, (iii) the broken line between 
d m + 1 points > ei — >• ei + e 2 — > ... — > ei + . . . + ed m (where ej is the 
i canonical basis vector in R dm , i.e., all of its coordinates are zero except the 
i, which is 1), and (iv) the skeleton of the unit square. Thus, the number of 
components M was equal to 4, and the dimension of the components (d m ) can 
be different and scaled. For illustration, see Fig. 1(b). In the ikeda test, hidden 
s™ = [s™!, s™ 2 ] €E K 2 sources realized the ikeda map 

s? +1>1 = 1 + A ro [<i cos«) - 8% sin«)], (9) 

ST + 1,2 = Amid Sin «) + C2 COS«)], (10) 

where A m is a parameter of the dynamical system and u>™ = 0.4— 1+ ^m )2+( s ™ y ■ 

M = 2 was chosen with initial points s{ = [20; 20], = [—100, 30] and parame- 
ters Ai = 0.9994, A 2 = 0.998, see Fig. 1(c) for illustration. 
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(a) (b) (c) 

Fig. 1: Illustration of the (a) smiley, (b) d-geom and (c) ikeda databases. 



4.2 Performance Measure, the Amari-index 

The identification of the fAR-IPA model is ambiguous, the hidden s m sources 
can be estimated up to ISA ambiguities. These ambiguities are however simple 
[12]: the components of equal dimension can be recovered up to the permutation 
(of equal dimension) and invertible transformation within the subspaces. Thus, 
in the ideal case, the product of the ISA demixing matrix Wisa and the ISA 



See http://www.smileyworld.com. 
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mixing matrix A, G = WisaA is a block-permutation matrix. This property 
can be measured by a simple extension of the Amari-index [13]. Namely, one can 
(i) assume without loss of generality that the component dimensions and their 
estimations are ordered in increasing order (di < ... < 4m, di < • ■ ■ < ^m), 



(ii) decompose G into di x dj blocks (G = [G u ] , 



and define g %3 as the 



sum of the absolute values of the elements of the matrix G' J € 



Then 



the Amari-index adapted to the ISA task of different component dimensions is 
defined as 



r(G) 



1 



2M(M - 1) 



E 



M 



M 

E 



maxj g %3 



(11) 



One can see that < r(G) < 1 for any matrix G, and r(G) = if and only if G 
is block-permutation matrix with di x dj sized blocks. r(G) = 1 is in the worst 
case, i.e, when all the g 13 elements are equal in absolute value. 



4.3 Simulations 

Results on databases smiley, d-geom and ikeda are provided here. For illustration 
purposes, we chose fAR order p = 1 and used the recursive Nadaraya- Watson 
(8) for functional AR estimation with the Gaussian kernel. The ISA subprob- 
lem was solved on the basis of the ISA separation theorem: the estimated ICA 
elements were clustered. The kernel canonical correlation technique [14] was ap- 
plied to estimate the dependence of the ICA elements. The permutation search 
(clustering step) was carried out by greedy optimization for tasks of known 
component dimensions (smiley, d-geom datasets). We employed the NCut [15] 
spectral technique on the ikeda dataset to estimate unknown dimensions and 
to perform clustering. FastICA [16] was used for the ICA estimation. Mixing 
matrix A was random orthogonal. For dataset smiley and d-geom, f was the 
composition of a random F matrix with entries distributed uniformly on in- 
terval [0, 1] and the noninvertible sine function. The Amari-index (Section 4.2) 
was used to evaluate the performance of the proposed fAR-IPA method. For 
each individual parameter, 10 random runs were averaged. Our parameters in- 
cluded T, the sample number of observations x 4 , and bandwith j3 € (0, 1/D) to 
study the robustness of the kernel regression approach. (3 was reparametcrized 
as /3 = and f3 c was chosen from the set ^, Jq, The perfor- 

mance of the method is summarized by notched boxed plots, which show the 
quartiles (Q\, Qi, Q%), depict the outliers, i.e., those that fall outside of interval 
[Qi — 1.5(Qs — Q\),Qz + 1.5(Q3 — Qi)] by circles, and whiskers represent the 
largest and smallest non-outlier data points. 

Figure 2 demonstrates that the algorithm was able to uncover the hidden 
components with high precision for the smiley dataset. Figure 2(a) illustrates 
the M — 2 (D = 4) case, Fig. 2(b) indicates that the problem with M = 6 
components (D = 12) for T = 50, 000 — 100, 000 samples is still amenable to the 
method. According to the figures, the estimation is robust with respect to the 
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choice of bandwith. The obtained source estimations are illustrated in Fig. 2(c)- 
(e). 
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Fig. 2: Illustration of the estimations on the smiley dataset. (a)-(b): Amari-index 
as a function of the sample number, for M = 2 and M = 6 components, re- 
spectively. The estimation error is plotted on log scale for different bandwith 
parameters, (c): observed signal x. (c): estimated components (e m ) with aver- 
age (closest to the median) Amari-index for M = 6, /3 C = T = 100, 000. (d): 
Hinton-diagram of matrix G for (e)-it is approximately a block-permutation 
matrix with 2x2 blocks. 



Our experiences concerning the d-geom and the ikeda datasets are summa- 
rized in Fig. 3. In accordance with the smiley test, the dimension of the d-geom 
problem was D — 12, however the dimensions of the hidden components were 
different and unknown to the algorithm: d\ = 2, d2 — d% = 3, d± = 4. As it can 
be seen from Fig. 3(a) the method provides precise estimations on the d-geom 
database for sample number T = 100, 000 — 150, 000. Hinton-diagram of matrix 
G with average (closest to the median) Amari-index is depicted in Fig. 3(c). 
Our third example is the ikeda database. As it is illustrated in Fig. 3(b), in this 
case an autoregressive approximation (AR-IPA) could not find the proper sub- 
spaces. Nevertheless, the Amari-index values of Fig. 3(b) show that a functional 
AR-IPA approach was able to recover the hidden subspaces, for sample number 
T > 10, 000. The figure also shows that the estimation is precise for a wide range 
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of bandwith parameters. Hidden sources with average Amari-indcx uncovered by 
the method are illustrated Fig. 3(d)-(f). 
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Fig. 3: Illustration of the estimations on the d-geom and the ikeda dataset. (a)- 
(b): Amari-index on log scale as a function of the sample number for different 
bandwith parameters on the d-geom (with component dimensions: d% — 2, di = 
d% — 3, d^ — 4) and the ikeda database, respectively, (c): Hinton-diagram of G 
with average (closest to the median) Amari-index for dataset d-geom, (3 C = 
T = 150,000-it is approximately a block-permutation matrix with one 2x2, 
two 3x3 and one 4x4 block, (d)-(f): estimation with average Amari-index 
for database ikeda, f3 c = |, T = 20,000. (d): observation x. (f): estimated 
components (s m ). (e): Hinton-diagram of matrix G for (f)-it is approximately 
a block-permutation matrix with 2x2 blocks. 



5 Conclusions 

In this paper we (i) extended independent subspace analysis (ISA) to the 
not strictly stationary domain, (ii) relaxed the constraint of decoupled (block- 
decorrelated) dynamics, and (iii) simultaneously addressed the case of unknown 
source component dimensions. The temporal evolution of the sources was cap- 
tured by functional autoregressive (fAR) processes. We generalized the ISA sepa- 
ration technique to the derived fAR setting (fAR-IPA, IPA- independent process 
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analysis) and reduced the solution of the problem to fAR identification and ISA. 
The fAR estimation was carried out by the Nadaraya- Watson kernel regres- 
sion method with strong consistency guarantee. We extended the Amari-indcx 
to different dimensional components and illustrated our technique by numeri- 
cal experiments. According to the experiences, the fAR-IPA identification can 
be accomplished robustly and can be advantageous compared to a parametric 
approach. The robustness of the separation principle indicate that it can be 
extended to a larger class of processes. 
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